# 0. initial setup ----
# clear workspace 
rm(list = ls()) 

# set seed
set.seed(11081984)

# set working directory
setwd(YOUR WORKING DIRECTORY)

# define background color
bg_colour <- "white"

# create folder to export figures and tables
dir.create("figures and tables")


## 0.1 define custom functions ----
main_model <- function(mydata, mymodel, myx_limit, myx_break) {
  ggplot(data = mydata %>%
           filter(model == mymodel)) + 
    geom_point(aes(estimate, factor), size = 4) +
    geom_point(data = mydata %>% filter(model == mymodel & estimate == 0), 
               aes(estimate, factor), size = 8, colour = "white") +
    geom_linerange(aes(y = factor, xmin = lower_95, xmax = upper_95), 
                   linewidth = 1.1) +
    geom_vline(xintercept = 0, 
               linetype = "dashed") +
    scale_x_continuous(name = "\nCoefficent estimates",
                       limits = c(-1 * myx_limit - 0.35,
                                  myx_limit + 0.35), 
                       breaks = c(-1 * myx_break, 
                                  0, 
                                  myx_break),
                       labels = c(paste("\n -", 
                                        as.character(myx_break), 
                                        "\n (decreases\n satisfaction)"),
                                  "0",
                                  paste("\n ", 
                                        as.character(myx_break), 
                                        "\n (increases\n satisfaction)"))) +
    scale_y_discrete(name = paste("Model", as.character(mymodel), "\n")) +
    theme_classic(base_size = 40) +
    theme(plot.background = element_rect(fill = bg_colour),
          axis.ticks.y=element_blank())
} 
check_install_load <- function() {
  installed_packs <- rownames(installed.packages())
  packs_to_install <- required_packs[!(required_packs %in% installed_packs)]
  if (length(packs_to_install) > 0) {
    install.packages(packs_to_install)
  }
  lapply(required_packs, require, character.only = TRUE)  
}
pos_scatter <- function(df, voters, issue_1, issue_2) {
  
  df <- df %>% filter(voter_party == voters)
  df_wdp <- df %>% filter(treatment == 3)
  party_abb <- case_when(voters == 1 ~ "PVV",
                         voters == 2 ~ "CDA",
                         voters == 3 ~ "D66",
                         voters == 4 ~ "GroenLinks",
                         voters == 5 ~ "SP",
                         voters == 6 ~ "PvdA",
                         voters == 7 ~ "CU",
                         voters == 8 ~ "PvdD")
  party_abb <- paste0(party_abb," ","(N=", dim(df)[1], ")")
  logroll1_x <- mean(df$ppi_1)
  logroll1_y <- mean(df$cppi_2)
  logroll2_x <- mean(df$cppi_1)
  logroll2_y <- mean(df$ppi_2)
  compromise_x <- mean(df_wdp$wdpi_1)
  compromise_y <- mean(df_wdp$wdpi_2)
  mean_voter_x <- mean(df$vpi_1)
  mean_voter_y <- mean(df$vpi_2)
  label_x <- case_when(issue_1 == 1 ~ "\nTax rate for the highest\nearners should be higher",
                       issue_1 == 2 ~ "\nRetirement age must rise\nas life expectancy increases",
                       issue_1 == 3 ~ "\nTo reduce nitrogen emissions,\nlivestock must be substantially reduced",
                       issue_1 == 4 ~ "\nSchools may no longer refuse\nstudents on the basis of religion",
                       issue_1 == 5 ~ "\nThe Netherlands must\nreceive more refugees",
                       issue_1 == 6 ~ "\nJudges should impose prison sentences\ninstead of community service more often",)
  label_y <- case_when(issue_2 == 1 ~ "Tax rate for the highest\nearners should be higher\n",
                       issue_2 == 2 ~ "Retirement age must rise\nas life expectancy increases\n",
                       issue_2 == 3 ~ "To reduce nitrogen emissions,\nlivestock must be substantially reduced\n",
                       issue_2 == 4 ~ "Schools may no longer refuse\nstudents on the basis of religion\n",
                       issue_2 == 5 ~ "The Netherlands must\nreceive more refugees\n",
                       issue_2 == 6 ~ "Judges should impose prison sentences\ninstead of community service more often\n",)
  
  fig <- ggplot(data = df) +
    geom_jitter(aes(vpi_1, vpi_2),
                width = 0.2, height = 0.2,
                shape = 1,
                size = 6) +
    geom_point(aes(x = logroll1_x, y = logroll1_y), 
               shape = 18,
               size = 9) +
    geom_point(aes(x = logroll2_x, y = logroll2_y), 
               shape = 18,
               size = 9) +
    geom_point(aes(x = compromise_x, y = compromise_y), 
               shape = 15,
               size = 9) +
    geom_vline(xintercept = mean_voter_x, 
               linetype = "dashed") +
    geom_hline(yintercept = mean_voter_y, 
               linetype = "dashed") +
    scale_x_continuous(name = label_x,
                       breaks = seq(1,5,1),
                       labels = c("strongly\nagree",
                                  "agree",
                                  "neutral",
                                  "disagree",
                                  "strongly\ndisagree")) +
    scale_y_continuous(name = label_y,
                       breaks = seq(1,5,1),
                       labels = c("strongly\nagree",
                                  "agree",
                                  "neutral",
                                  "disagree",
                                  "strongly\ndisagree")) +
    ggtitle(party_abb) +
    theme_classic(base_size = 30) +
    theme(plot.background = element_rect(fill = "white"))
  
  return(fig)
  
}


## 0.2 define and load required packages
required_packs <- c("tidyverse",
                    "stargazer",
                    "modelsummary",
                    "mediation",
                    "ggeffects",
                    "glmnet",
                    "margins",
                    "multcomp")
check_install_load()



# 1. empirical analysis ----
## 1.1 load analysis data set ----
load("data/analysis_data.Rdata")


## 1.2 restrict analysis to non-VVD supporters (see footnote 12) ----
df_analysis_full <- df_analysis
df_analysis <- df_analysis %>%
  filter(voter_party != 9) %>%
  mutate(deal_logroll = as.factor(deal_logroll))


## 1.3 define DV, IVs and control variables ----
# outcome variable
outcome <- "Q28_voter_deal_satis"

# independent variables 
ind_var1 <- "deal_logroll"
ind_var2 <- "deal_logroll * voter_deal_wdis"

# control variables
controls1 <- c("voter_deal_wdis",
              "hostility",
              "as.factor(pid_strength)",
              "as.factor(party_dyad)")
controls2 <- c("hostility",
               "as.factor(pid_strength)",
               "as.factor(party_dyad)")


## 1.4 estimate models for satisfaction with government deal ----
# model 1 - hypothesis 1
model_1 <- lm(paste(outcome, paste(c(ind_var1, controls1), 
                                   collapse = "+"), 
                    sep = "~"),
              data = df_analysis)

# model 2 - hypothesis 2
model_2 <- lm(paste(outcome, paste(c(ind_var2, controls2), 
                                   collapse = "+"), 
                    sep = "~"),
              data = df_analysis)



# 2. preregistered empirical analyses ----
## 2.1 figure 4 - determinants of satisfaction with government deal ----
# extract coefficients from models 1 and 2
df_coefs_mod1 <- broom::tidy(model_1) %>% 
  mutate(model = 1)
df_coefs_mod2 <- broom::tidy(model_2) %>% 
  mutate(model = 2)
df_coefs_mod <- rbind(df_coefs_mod1, df_coefs_mod2) %>% 
  filter(term %in% c("deal_logrollTRUE",
                     "voter_deal_wdis",
                     "hostility",
                     "as.factor(pid_strength)2",
                     "as.factor(pid_strength)3",
                     "as.factor(pid_strength)4",
                     "deal_logrollTRUE:voter_deal_wdis")) %>%
  mutate(factor = recode(term,
                         "hostility" = "hostility",
                         "voter_deal_wdis" = "policy\ndistance\n(\u03B1)",
                         "deal_logrollTRUE" = "logroll\n(\u03B2)",
                         "deal_logrollTRUE:voter_deal_wdis" = "logroll \u00d7 policy\ndistance\n(\u03B3)",
                         "as.factor(pid_strength)2" = "convinced\nadherent",
                         "as.factor(pid_strength)3" = "not convinced\nadherent",
                         "as.factor(pid_strength)4" = "not adherent")) %>%
  dplyr::select(factor, estimate, std.error, model) %>% 
  mutate(lower_95 = estimate - 1.96 * std.error,
         upper_95 = estimate + 1.96 * std.error) %>%
  add_row(factor = "very convinced\nadherent (ref.)",
          estimate = 0,
          model = 1,
          .before = 4) %>%
  add_row(factor = "very convinced\nadherent (ref.)",
          estimate = 0,
          model = 2,
          .before = 4) %>%
  add_row(factor = " ",
          estimate = 0,
          model = 1,
          .before = 4) %>%
  mutate(factor = factor(factor, levels = c("not adherent",
                                            "not convinced\nadherent",
                                            "convinced\nadherent",
                                            "very convinced\nadherent (ref.)",
                                            "hostility",
                                            " ",
                                            "logroll \u00d7 policy\ndistance\n(\u03B3)",
                                            "logroll\n(\u03B2)",
                                            "policy\ndistance\n(\u03B1)")))
rm(df_coefs_mod1, df_coefs_mod2)

# define limits and breaks for x-axis based on coefficient estimates  
x_limit <- round(max(abs(df_coefs_mod[,4:5]), na.rm = TRUE),1)
x_break <- round(max(abs(df_coefs_mod[,4:5]), na.rm = TRUE))

# plot left and right panel of figure 3
fig_mod1 <- main_model(df_coefs_mod,
                       1,
                       x_limit,
                       x_break)

fig_mod2 <- main_model(df_coefs_mod,
                       2,
                       x_limit,
                       x_break)

# export left and right panel from figure 4 
jpeg("figures and tables/fig_coef_model1.jpeg", quality = 100,
     height = 1200, 
     width = 1000,
     units = "px")
fig_mod1
dev.off()

jpeg("figures and tables/fig_coef_model2.jpeg", quality = 100, 
     height = 1200, 
     width = 1000,
     units = "px")
fig_mod2
dev.off()


## 2.2 figure 5 - satisfaction with government deal conditional on policy distance ----
# generate predicted values based on model 2
df_pred_mod2 <- ggpredict(model_2, terms = c("voter_deal_wdis",
                                             "deal_logroll")) %>%
  as.data.frame() %>%
  mutate(group = recode(group, 
                        "FALSE" = "compromise", 
                        "TRUE" = "logroll"),
         x = case_when(group == "compromise" ~ x - 0.025,
                       group == "logroll" ~ x + 0.025))

# marginal effect of a log-roll conditional on policy distance
voter_deal_wdis <- seq(0,4,0.1)
mat_logroll <- matrix(c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 1)
mat_test_H2 <- mat_logroll %x% rep(1, length(voter_deal_wdis))
mat_test_H2[,15] <- voter_deal_wdis

# estimate marginal log-roll effect 
test_H2 <- multcomp::glht(model_2, mat_test_H2)
df_marg_mod2 <- summary(test_H2)$test
df_marg_mod2 <- as.data.frame(cbind(df_marg_mod2$coefficients, 
                                    df_marg_mod2$sigma, voter_deal_wdis))
colnames(df_marg_mod2) <- c("estimate", 
                            "std.error", 
                            "values_dis")
df_marg_mod2 <- df_marg_mod2 %>% 
  mutate(sig = case_when(estimate/std.error < qnorm(p = .05/2, lower.tail = TRUE) ~ 1,
                         estimate/std.error > qnorm(p = .05/2, lower.tail = FALSE) ~ 2,
                         abs(estimate) <= abs(std.error) * qnorm(p = .05/2, lower.tail = FALSE) ~ 0)) %>%
  filter(sig == 0) %>%
  summarise(min_dis = min(values_dis),
            max_dis = max(values_dis))

# share of respondents with significant positive/negative log-roll effect
pos <- round(table(df_analysis$voter_deal_wdis < df_marg_mod2$min_dis)[2]/nrow(df_analysis)*100,1)
neg <- round(table(df_analysis$voter_deal_wdis > df_marg_mod2$max_dis)[2]/nrow(df_analysis)*100,1)

# plot predicted values over policy distance
fig_pred_mod2 <- ggplot(df_pred_mod2, aes(x, predicted)) +
  geom_point(aes(shape = group), size = 4) +
  scale_shape_manual(values = c(19, 0)) +
  geom_linerange(aes(x = x, ymin = conf.low, ymax = conf.high), 
                 linewidth = 1.1) +
  geom_vline(xintercept = df_marg_mod2$min_dis) + 
  geom_vline(xintercept = df_marg_mod2$max_dis) +
  scale_y_continuous(name = "Satisfaction with government deal\n",
                     limits = c(0,10), 
                     breaks = c(0,5,10),
                     labels = c("0 (break off\nnegotiations)",
                                "5",
                                "10 (sign the\nagreement)")) +
  scale_x_continuous(name = "\nPolicy distance",
                     limits = c(-0.1,4.1),
                     breaks = seq(0,4,1),
                     labels = c("0\n(close to\nideal point)",
                                "1", "2", "3",
                                "4\n(far from\nideal point)")) +
  geom_segment(aes(x = df_marg_mod2$min_dis - 0.05, xend = 0 + 0.05,
                   y = 1, yend = 1),
               linewidth = 1.1,
               arrow = arrow(length = unit(0.5, "cm"))) +
  geom_segment(aes(x = df_marg_mod2$max_dis + 0.05, xend = 4 - 0.05,
                   y = 1, yend = 1),
               linewidth = 1.1,
               arrow = arrow(length = unit(0.5, "cm"))) + 
  geom_text(x = 0.5, 
            y = 1.8, 
            label = paste0("significant positive\nlogroll effect (",
                           pos,
                           " %)"),
            size = 8,
            colour = "black") + 
  geom_text(x = 3, 
            y = 1.8, 
            label = paste0("significant negative\nlogroll effect (",
                           neg,
                           " %)"),
            size = 8,
            colour = "black") + 
  theme_classic(base_size = 40) +
  theme(plot.background = element_rect(fill = bg_colour),
        legend.title = element_blank(),
        legend.background = element_rect(fill = bg_colour))

# export figure 5
jpeg("figures and tables/fig_pred_model2.jpeg", quality = 100, 
     height = 1000, 
     width = 1600,
     units = "px")
fig_pred_mod2
dev.off()



# 3. additional exploratory analyses ----
## 3.1 coherent vs incoherent log-rolls ----
df_analysis <- df_analysis %>%
  mutate(coherent_logroll = case_when((ppi_1 == 1 | ppi_1 == 5) & treatment == 1 ~ 1,
                                      (ppi_2 == 1 | ppi_2 == 5) & treatment == 2 ~ 1,
                                      treatment == 3 ~ 0),
         coherent_logroll = replace_na(coherent_logroll, -1),
         # deviate from empirical coding rule to accommodate CDA-VVD dyad
         coherent_logroll = if_else(voter_party == 2 & coherent_logroll == -1,
                                    1, coherent_logroll),
         coherent_logroll = relevel(as.factor(coherent_logroll), ref = 2))
table(df_analysis$deal_logroll, 
      df_analysis$coherent_logroll)

# estimate model 3
model_3 <- lm(paste(outcome, 
                    paste(c("as.factor(coherent_logroll)", 
                            controls1), 
                          collapse = "+"), sep = "~"),
              data = df_analysis) 


## 3.2 figure 6 - empirical effects of coherent vs incoherent log-rolls ----
# extract coefficients of model 3
df_coefs_mod3 <- broom::tidy(model_3) %>% 
  filter(term %in% c("as.factor(coherent_logroll)-1",
                     "as.factor(coherent_logroll)1",
                     "voter_deal_wdis",
                     "hostility",
                     "as.factor(pid_strength)2",
                     "as.factor(pid_strength)3",
                     "as.factor(pid_strength)4")) %>%
  mutate(factor = recode(term,
                         "hostility" = "hostility",
                         "voter_deal_wdis" = "policy\ndistance",
                         "as.factor(pid_strength)2" = "convinced\nadherent",
                         "as.factor(pid_strength)3" = "not convinced\nadherent",
                         "as.factor(pid_strength)4" = "not adherent",
                         "as.factor(coherent_logroll)-1" = "incoherent\nlogroll",
                         "as.factor(coherent_logroll)1"  = "coherent\nlogroll")) %>%
  dplyr::select(factor, estimate, std.error) %>% 
  mutate(lower_95 = estimate - 1.96 * std.error,
         upper_95 = estimate + 1.96 * std.error) %>%
  add_row(factor = "very convinced\nadherent (ref.)",
          estimate = 0, 
          .before = 5) %>%
  add_row(factor = "compromise\n(ref.)",
          estimate = 0, 
          .before = 7) %>%
  mutate(factor = factor(factor, levels = c("not adherent",
                                            "not convinced\nadherent",
                                            "convinced\nadherent",
                                            "very convinced\nadherent (ref.)",
                                            "hostility",
                                            "coherent\nlogroll",
                                            "compromise\n(ref.)",
                                            "incoherent\nlogroll",
                                            "policy\ndistance")),
         model = 3)

# plot figure coefficient estimates of model 3
fig_mod3 <- main_model(df_coefs_mod3,
                       3,
                       x_limit,
                       x_break)

# export figure 6
jpeg("figures and tables/fig_coef_model3.jpeg", quality = 100, 
     height = 1000, 
     width = 1600,
     units = "px")
fig_mod3
dev.off()


## 3.3 mediation analysis - would you vote for same party in the next election? ----
df_analysis_vote <- df_analysis %>% 
  filter(!is.na(Q29_vote_again)) %>%
  mutate(logroll = if_else(deal_logroll == TRUE, 1 ,0))

model_mediator <- lm(Q28_voter_deal_satis ~ 
                       logroll * voter_deal_wdis +
                       hostility +
                       as.factor(pid_strength) +
                       as.factor(party_dyad),
                     data = df_analysis_vote)

model_outcome <- glm(Q29_vote_again ~ 
                       Q28_voter_deal_satis * voter_deal_wdis + 
                       logroll * voter_deal_wdis +
                       hostility +
                       as.factor(pid_strength) +
                       as.factor(party_dyad),
                     data = df_analysis_vote,
                     family = binomial)

model_mediation <- mediate(model_mediator, 
                           model_outcome,
                           treat = "logroll",
                           mediator = "Q28_voter_deal_satis",
                           sims = 1000,
                           covariates = list(voter_deal_wdis = 0))
summary(model_mediation)

model_mod_med <- test.modmed(model_mediation, 
                             covariates.1 = list(voter_deal_wdis = 0),
                             covariates.2 = list(voter_deal_wdis = 4))

model_mod_med


## 3.4 figure 7 - effect of dissatisfaction with government deal on probability to vote for same party ----
# generate predicted probabilities
df_coefs_mod4 <- ggpredict(model_outcome, 
                           terms = "Q28_voter_deal_satis",
                           condition = c(logroll = 1,
                                         voter_deal_wdis = median(df_analysis_vote$voter_deal_wdis),
                                         hostility = median(df_analysis_vote$hostility),
                                         pid_strength = 2,
                                         party_dyad = 1))
margins_summary(model_outcome, 
                variables = "Q28_voter_deal_satis")

# plot predicted probability to vote for same party again
fig_pred_mod4 <- ggplot(df_coefs_mod4, aes(x, predicted)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1) +
  scale_y_continuous(name = "Predicted probability to \n vote for same party again \n",
                     limits = c(0,1),
                     breaks = c(0, 0.5, 1)) +
  scale_x_continuous(name = "\n Satisfaction with government deal",
                     limits = c(0,11), 
                     breaks = c(0,5,10),
                     labels = c("0\n(break off\nnegotiations)",
                                "5",
                                "10\n(sign the\nagreement)")) +
  geom_rug(data = df_analysis_vote, aes(voter_deal_satis_jitter, Q29_vote_again),
           alpha = .8, position = "identity", sides = "b") + 
  theme_classic(base_size = 40) +
  theme(plot.background = element_rect(fill = bg_colour))

# export figure 7
jpeg("figures and tables/fig_pred_model4.jpeg", quality = 100, 
     height = 1000, 
     width = 1600, 
     units = "px")
fig_pred_mod4
dev.off()



# 4. appendix A - descriptive statistics and empirical results ----
## 4.1 table A2 - descriptive statistics ----
df_descriptive <- df_analysis %>% 
  dplyr::select(Q28_voter_deal_satis,
                voter_deal_wdis,
                hostility,
                treatment,
                pid_strength) %>% 
  drop_na() %>%
  as.data.frame()

round(table(df_descriptive$treatment)/sum(table(df_descriptive$treatment))*100, 2)
round(table(df_descriptive$pid_strength)/sum(table(df_descriptive$pid_strength))*100, 2)

stargazer(df_descriptive, 
          nobs = FALSE,
          median = TRUE,
          digits = 2,
          type = "html", 
          out= "figures and tables/table_descriptives.doc")


## 4.2 figure A1 - voter policy positions ----
figa1 <- list()  
figa1[[1]] <- pos_scatter(df_analysis, 
                          1, 1, 2)
figa1[[2]] <- pos_scatter(df_analysis, 
                          2, 1, 4)
figa1[[3]] <- pos_scatter(df_analysis, 
                          3, 3, 6)
figa1[[4]] <- pos_scatter(df_analysis, 
                          4, 3, 6)
figa1[[5]] <- pos_scatter(df_analysis, 
                          5, 5, 2)
figa1[[6]] <- pos_scatter(df_analysis, 
                          6, 1, 6)
figa1[[7]] <- pos_scatter(df_analysis, 
                          7, 4, 5)
figa1[[8]] <- pos_scatter(df_analysis, 
                          8, 3, 6)

for(i in 1:8) {
  name <- paste("figures and tables/fig_voters_party_", i, ".jpeg", sep = "")
  
  jpeg(name, quality = 100, 
       height = 1000, 
       width = 1000,
       units = "px")
  print(figa1[[i]])
  dev.off()
  
}


## 4.3 table A3 - model results ----
stargazer(model_1, model_2,
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          omit.stat = "f",
          type = "html",
          out = "figures and tables/table_model_results.doc")


## 4.4 table A4 - exploratory analyses ----
stargazer(model_3, model_outcome, 
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          omit.stat = "f",
          type = "html",
          out = "figures and tables/table_exploratory_results.doc")



# 5. appendix B - robustness of empirical findings ----
## 5.1 models R1.X - excluding respondents post 25/3 ----
df_analysis_sub <- df_analysis %>% 
  filter(id %in% analysis_sub)

# share of respondents, see footnote 5
nrow(df_analysis_sub)/nrow(df_analysis)

mR1.1 <- lm(paste(outcome, paste(c(ind_var1, controls1), 
                                  collapse = "+"), sep = "~"),
             data = df_analysis_sub)

mR1.2 <- lm(paste(outcome, paste(c(ind_var2, controls2), 
                                  collapse = "+"), sep = "~"),
             data = df_analysis_sub)


## 5.2 models R2.X - using binary policy positions ----
ind_var1_R2 <- ind_var1
ind_var2_R2 <- "deal_logroll * voter_deal_dis_bin"

controls1_R2 <- c("voter_deal_dis_bin",
                  "hostility",
                  "as.factor(pid_strength)",
                  "as.factor(party_dyad)")
controls2_R2 <- controls2

mR2.1 <- lm(paste(outcome, paste(c(ind_var1_R2, controls1_R2), 
                                  collapse = "+"), sep = "~"),
             data = df_analysis)

mR2.2 <- lm(paste(outcome, paste(c(ind_var2_R2, controls2_R2), 
                                  collapse = "+"), sep = "~"),
             data = df_analysis)


## 5.3 models R3.X - using non-weighted average positions ----
ind_var1_R3 <- ind_var1
ind_var2_R3 <- "deal_logroll * voter_deal_dis"

controls1_R3 <- c("voter_deal_dis",
                  "hostility",
                  "as.factor(pid_strength)",
                  "as.factor(party_dyad)")
controls2_R3 <- controls2

mR3.1 <- lm(paste(outcome, paste(c(ind_var1_R3, controls1_R3), 
                                  collapse="+"), sep="~"),
             data = df_analysis)

mR3.2 <- lm(paste(outcome, paste(c(ind_var2_R3, controls2_R3), 
                                  collapse="+"), sep="~"),
             data = df_analysis)


## 5.4 models R4.X - using Euclidean policy distance ----
ind_var1_R4 <- ind_var1
ind_var2_R4 <- "deal_logroll * voter_deal_edis"

controls1_R4 <- c("voter_deal_edis",
                  "hostility",
                  "as.factor(pid_strength)",
                  "as.factor(party_dyad)")
controls2_R4 <- controls2

mR4.1 <- lm(paste(outcome, paste(c(ind_var1_R4, controls1_R4), 
                                 collapse = "+"), sep = "~"),
            data = df_analysis)

mR4.2 <- lm(paste(outcome, paste(c(ind_var2_R4, controls2_R4), 
                                 collapse = "+"), sep = "~"),
            data = df_analysis)


## 5.5 models R5.X - using full sample including VVD supporters ----
mR5.1 <- lm(paste(outcome, paste(c(ind_var1, controls1), 
                                  collapse = "+"), sep = "~"),
                   data = df_analysis_full)

mR5.2 <- lm(paste(outcome, paste(c(ind_var2, controls2), 
                                  collapse = "+"), sep = "~"),
                   data = df_analysis_full)


## 5.6 models R6.X - exploring null-effect of policy distance in the compromise condition ---- 
mR6.1 <- glm(Q29_vote_again ~ 
               Q28_voter_deal_satis * logroll + 
               voter_deal_wdis +
               hostility +
               as.factor(pid_strength) +
               as.factor(party_dyad),
             data = df_analysis_vote,
             family = binomial)

margins_summary(mR6.1, 
                variables = "Q28_voter_deal_satis",
                at = list(logroll = c(0, 1)))

df_analysis <- mutate(df_analysis,
                      voter_pos_centr = 
                        (((abs(Q11_vpi_taxes - 3) + 
                           abs(Q11_vpi_pensions - 3) +
                           abs(Q11_vpi_environ - 3) +
                           abs(Q11_vpi_education - 3) +
                           abs(Q11_vpi_immigration - 3) +
                           abs(Q11_vpi_interior - 3)) * -1)/12) + 1,
                      voter_pos_extr = (vpi_1 - 3)^2 + (vpi_2 - 3)^2,
                      voter_centrist = if_else(voter_pos_extr <= 1, "centrist", "extreme"))

cor.test(df_analysis$voter_pos_centr, df_analysis$Q5_pol_interest)
cor.test(df_analysis$voter_pos_centr, df_analysis$Q6_campaign)
cor.test(df_analysis$Q5_pol_interest, df_analysis$Q6_campaign)


## 5.7 figure B.1 ----
df_analysis_cen <- filter(df_analysis,
                          voter_deal_wdis < quantile(voter_deal_wdis, 0.875) &
                            voter_deal_wdis > quantile(voter_deal_wdis, 0.05)) %>%
  mutate(deal_logroll = if_else(deal_logroll == TRUE, "logroll", "compromise"))

fig_b1 <- ggplot(df_analysis_cen,
       aes(y = Q28_voter_deal_satis,
           x = voter_deal_wdis)) +
  geom_jitter(data = df_analysis_cen,
              aes(color = deal_logroll),
              width = 0.05, 
              height = 0.5,
              size = 4) +
  geom_smooth(aes(linetype = deal_logroll),
              method = lm, formula = y ~ splines::bs(x, 7), se = FALSE,
              linewidth = 1.5) +
  scale_color_manual(values = c(compromise = "grey",
                                logroll = "black")) + 
  scale_linetype_manual(values = c(compromise = 2,
                                   logroll = 1)) +
  scale_y_continuous(name = "Satisfaction with government deal\n",
                     limits = c(0,10), 
                     breaks = c(0,5,10),
                     labels = c("0 (break off\nnegotiations)",
                                "5",
                                "10 (sign the\nagreement)")) +
  scale_x_continuous(name = "\nPolicy distance") +
  theme_classic(base_size = 40) +
  theme(plot.background = element_rect(fill = bg_colour),
        legend.key.size = unit(5,"line"),
        legend.title=element_blank())

jpeg("figures and tables/fig_b1.jpeg", quality = 100, 
     height = 1000, 
     width = 1600, 
     units = "px")
fig_b1
dev.off()


## 5.8 table B1 - alternative model specifications ----
stargazer(mR1.1, mR1.2,
          mR2.1, mR2.2,
          mR3.1, mR3.2,
          mR4.1, mR4.2,
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          omit.stat = "f",
          type = "html",
          out = "figures and tables/table_robustness.doc")

stargazer(mR5.1, mR5.2,
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          omit.stat = "f",
          type = "html",
          out = "figures and tables/table_robustness_VVD.doc")